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We report both experimental measurements and molecular simulations of 
the melting and freezing behavior of simple fluids in porous media. The 
experimental studies are for carbon tetrachloride and nitrobenzene in 
controlled pore glass (CPG), Vycor and activated carbon fibers (ACF). 
Differential scanning calorimetry (DSC) was used to determine the melt- 
ing point in these porous materials. In the case of nitrobenzene (which has 
a dipole moment), dielectric spectroscopy was also used to determine 
melting points. Measurements by the two methods were in excellent 
agreement. The melting point was found to be depressed relative to the 
bulk value for both fluids in the silica based pores. While nitrobenzene in 
ACF showed a very small shift in melting temperature, a large increase in 
the melting temperatures was found for CCly in ACF's, relative to the 
bulk. With the exception of the smallest pores, the melting point depres- 
sion was proportional to the reciprocal of the pore diameter, in agreement 
with the Gibbs-Thomson equation. Structural information about the dif- 
ferent confined phases was obtained by measuring the dielectric relaxation 
times using dielectric relaxation spectroscopy. 

Monte Carlo simulations were used to determine the shift in the melting 
point, 7,, for simple fluids in pores having repulsive, weakly attractive 
and strongly attractive walls chosen to model a range of porous materials 
from silica based pores to carbon based pores, in terms of energy of the 
fluid-wall interaction. The strength of the interaction of the fluid with the 
pore wall is shown to have a large effect on the shift in T,,, with 7, being 
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reduced for weakly attracting walls. For strongly attracting walls, such as 
graphitic carbon, the melting point increases for slit-shaped pores. These 
trends are compared with the experimental trends for different fluids con- 
fined in different porous materials based on corresponding states theory. 
A new method for calculating the free energies of solids in pores based on 
Landau theory is used, and it is shown that the solid-liquid transition re- 
mains first order in these confined systems. 


1. INTRODUCTION 


Recently there has been growing interest in the study of fluid-solid transi- 
tions in porous materials. Important questions are the nature of the phase tran- 
sition (first order vs. continuous), the direction of shift in the melting tempera- 
ture T, nature and origin of hysteresis, dimensionality cross-over due to in- 
creasing confinement, structural changes of the condensed phases in the re- 
stricted pore geometries, effect on latent heats, etc. Improved understanding of 
confinement effects on freezing are essential in areas relating to lubrication, 
adhesion, fabrication of nanomaterials and nanotribology. In addition, these 
studies can provide insight into mechanisms involved in frost heaving and dis- 
tribution of pollutants in soil. Freezing in porous media has also been widely 
employed in the characterization of porous materials using the method of ther- 
moporometry. Some other possible applications are in the development of novel 
materials for energy storage, high speed vehicles and aircraft braking sys- 
tems [1,2]. 

Experiments on freezing that have used porous silica glass as the confine- 
ment medium have always resulted in a decrease in the freezing temperature, 
T,, as compared to the bulk. Freezing of oxygen in sol-gel glasses was studied 
by Warnock et. al. [3] by a sub-pico second optical technique. In this method, 
birefringence in optical pump pulses caused by the rotational motion of the 
molecules in the liquid was used to measure the subsequent molecular orienta- 
tional relaxation time. A change in the value of the relaxation time provided an 
indication of the freezing temperature. The freezing temperature in the confined 
system was always depressed as compared to the bulk; the shift was larger for 
smaller pores, and as large as 10 K for the smallest (20 nm) pore. Unruh and 
co-workers [4] examined the melting behavior of indium metal in porous silica 
glasses by differential scanning calorimetry (DSC) measurements, and reported 
a large depression in melting point due to confinement. In view of the large 
body of experimental evidence for a decrease in the freezing temperature due to 
confinement ([5] and the refs. therein), it is tempting to assume that a decrease 
always occurs. However, a classical thermodynamic argument based on simple 
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capillary theory determines the freezing temperature as the point at which the 
chemical potential of the solid core inside the pore equals that of the surround- 
ing fluid (for example, see ref. [5]), 


AT, aS ay (Ys en Vwi )v 


re HA, 


(1) 


where T, is the bulk freezing temperature, y,, and y,, are the corresponding 
wall-solid and wall-fluid surface tensions, v is the molar volume of the liquid, 
Ap is the latent heat of melting in the bulk and H is the pore width. Equation 
(1) is sometimes referred to as the Gibbs-Thomson equation, and predicts that 
the sign of AT, = T, —T,, depends on whether y,,, is greater or less than y,,,. 
In a subsequent molecular simulation study of the effect of confinement on 
freezing of simple fluids in slit pores by Miyahara and Gubbins [6], it was 
shown that 7, was strongly affected by the strength of the attractive forces 
between the fluid molecules and the pore walls. For repulsive or weakly attrac- 
tive potentials, the shift in the freezing temperature AT, was negative. For 
strongly attracting walls such as carbons, an increase in T, was observed. 
Thus, Miyahara and Gubbins explained the disparate experimental trends on 
the direction of the shift in the freezing temperature and provided the connec- 
tion to the Gibbs-Thomson equation. The predictions of Miyahara and Gubbins 
were confirmed by free energy studies that provided the thermodynamic freez- 
ing temperature in confined systems; these studies established that the freezing 
transition remained first order despite the varied dimensionality [7,8]. 

In this article, we report a systematic study of the effect of the strength of 
the pore wall interaction on the phase behavior and molecular structure associ- 
ated with melting and freezing. The DSC and dielectric spectroscopy experi- 
ments involving porous glasses and ACF's were performed at Adama 
Mickiewicza University, Poznan, Poland, and Chiba University, Chiba, Japan 
and the molecular simulations were performed at North Carolina State Univer- 
sity, USA. 


2. METHODS 


Differential Scanning Calorimetry. Perkin-Elmer and Dupont thermal ana- 
lyzers were used to determine the melting temperatures and latent heats of fu- 
sion, by measuring the heat released in the melting of nitrobenzene and CCl, in 
the bulk as well as in the porous silica and activated carbon fibers. Temperature 
scanning rates of 0.5-2 K/min were used in the experiments. The background of 
each raw DSC spectrum was subtracted from the signal, based on a second- 
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order polynomial fit to the measured heat flow. The melting temperatures were 
determined from the position of the peaks of the heat flow signals and the latent 
heats were determined based on the scaled area under these signals. The melt- 
ing temperatures were reproducible to within 0.5°C for larger pores (225 nm); 
uncertainties were larger for the smaller pores. These uncertainties are a result 
of the width of the DSC peaks, which derive in part from variations in pore 
size, and geometry, and from the existence of metastable states. The latent heats 
were reproducible to within an accuracy of 5%. 


Dielectric relaxation spectroscopy. The relative permittivity of a medium, 
K` =k'-ik"', is in general a complex quantity whose real part x' (also known 
as the dielectric constant) is associated with the increase in capacitance due to 
the introduction of the dielectric. The imaginary component k"' is associated 
with mechanisms that contribute to energy dissipation in the system, due to 
viscous damping of the rotational motion of the dipolar molecules in alternating 
fields. The latter effect is frequency dependent. The experimental setup con- 
sisted of a parallel plate capacitor of empty capacitance C, = 4.2 pF. The ca- 
pacitance, C, and the tangent loss, tan(5), of the capacitor filled with nitro- 
benzene between the plates were measured using a Solartron 1260 gain imped- 
ance analyzer, in the frequency range 1 Hz - 10 MHz, for various temperatures. 
For the case of nitrobenzene in porous silica, the sample was introduced be- 
tween the capacitor plates as a suspension of 200 um mesh porous silica parti- 
cles in pure nitrobenzene. The contributions to x” were determined from, 
K'= hom K'= U3) 
ò K 


(2) 


In equation (2), C is the capacitance, C, is the capacitance without the dielec- 
tric and 6 is the angle by which current leads the voltage. Nitrobenzene was 
confined in porous silica (CPG and VYCOR), of pore widths H =50 nm to 
4nm at | atm. pressure. The freezing temperature in the bulk is 5.6°C (the liq- 
uid freezes to a monoclinic crystal). x° =K'-ix"', the complex dielectric per- 
mittivity, is measured as a function of temperature and frequency. 

For an isolated dipole rotating under an oscillating electric field in a viscous 
medium, the Debye dispersion relation is derived using classical mechanics, 


J 
KHSR 
1+i@t 


K =K' + (3) 


In equation (3), @ is the frequency of the applied potential and T is the orien- 
tational (rotational) relaxation time of a dipolar molecule. The subscript s re- 
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fers to static permittivity (low frequency limit, when the dipoles have sufficient 
time to be in phase with the applied field). The subscript œ refers to the optical 
permittivity (high frequency limit) and is a measure of the induced component 
of the permittivity. The dielectric constant is a natural choice of order parame- 
ter to study freezing of dipolar liquids, because of the large change in the ori- 
entational polarizability between the liquid and solid phases. The dielectric 
relaxation time was calculated by fitting the dispersion spectrum of the com- 
plex permittivity near resonance to the Debye model of orientational relaxation. 


Simulation. We performed Grand Canonical Monte Carlo (GCMC) simula- 
tions of Lennard-Jones fluids adsorbed in regular slit shaped pores of different 
pore widths H and varying fluid-wall strengths. Here H is the distance sepa- 
rating the planes through the centers of the surface-layer atoms on opposing 
pore walls. The interaction between the adsorbed fluid molecules is modeled 
using the Lennard-Jones (12,6) potential with size and energy parameters cho- 
sen to describe methane and CCl, [8,9]. The fluid-wall interaction is modeled 
using a ”10-4-3” Steele potential [10], 


E eE E KT A _(__% _ 
#.)=210.£,02.4 2 a ) ( > ) akal (4) 


Here, the o's and e's are the size and energy parameters in the LJ potential, 
the subscripts f and w denote fluid and wall respectively, z is the coordinate 
perpendicular to the pore walls and k, is the Boltzmann's constant. The fluid- 
wall interaction energy parameters corresponding to a graphite pore were taken 
from Ref. [10]. For a given pore width H , the total potential energy from both 
walls is given by, 


P pore (Z) = Pp (Z) + Pp (H-2) (5) 


The strength of the fluid wall interaction is determined by the parameter 
Q = 2mp,,€,,05,A. For example, in the case of methane, a = 0 corresponds to a 
purely repulsive or hard wall, @=0.76€, corresponds to a weakly attractive 
wall and @ = 2.0€,, corresponds to a strongly attractive graphite wall. 

The simulation runs were performed in the grand canonical ensemble, fix- 
ing the chemical potential u, the volume V of the pore and the temperature 
T . The system typically consisted of 600-700 adsorbed molecules. For the case 
of attractive pore-wall interaction, the adsorbed molecules formed distinct lay- 
ers parallel to the plane of the pore walls. A rectilinear simulation cell of 100, 
by 100, in the plane parallel to the pore walls was used, consistent with a 
cutoff of 5o, for the fluid-fluid interaction. The simulation was set up such 
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that insertion, deletion and displacement moves were attempted with equal 
probability, and the displacement step was adjusted to have a 50% probability 
of acceptance. Thermodynamic properties were averaged over 100-500 million 
individual Monte Carlo steps. The length of the simulation was adjusted such 
that a minimum of fifty times the average number of particles in the system 
would be inserted and deleted during a single simulation run. 


Free energy method. The method relies on the calculation of the Landau free 
energy as a function of an effective bond orientational order parameter ®, 
using GCMC simulations [8]. The Landau free energy A is defined by, 


A[{®]=—-k,7 In(P[®]) + constant (6) 


where P[®] is the probability of observing the system having an order pa- 
rameter value between ® and ®+6®. The probability distribution function 
P[®] is calculated in a GCMC simulation as a histogram, with the help of 
umbrella sampling. The grand free energy Q is then related to the Landau free 
energy by 


exp(—BQ) = | d@ exp(—BAL®]) (1) 


The grand free energy at a particular temperature can be calculated by numeri- 
cally integrating equation (7) over the order parameter space. We use a two- 
dimensional order parameter to characterize the order in each of the molecular 
layers. 


= > exp iss, ) 


b k=l 


p = 


Jj 


2 |(expti60, )), (8) 


P, measures the hexagonal bond order within each layer j. Each nearest 
neighbor bond has a particular orientation in the plane of the given layer, and is 
described by the polar coordinate 6. The index k runs over the total number 
of nearest neighbor bonds N, in layer j. The overall order parameter ® is an 
average of the hexagonal order in all the layers. We expect ®, = 0 when layer 
j has the structure of a two-dimensional liquid, ®, =1 in the crystal phase 
and 0< ®, <1 in an orientationally ordered layer. 


3. RESULTS 
3.1. Experiments and Simulations on Porous Silica 


DSC & Dielectric Spectroscopy. DSC scans corresponding to melting of ni- 
trobenzene are shown in Figure la. The deep minimum at 5. 6°C that is inde- 
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pendent of pore size corresponds to the melting of the bulk nitrobenzene in 
which the porous particles are suspended. In addition, a second minimum that is 
dependent on the pore size is also observed. This signal corresponds to the 
melting of nitrobenzene in the pores, and shifts to lower temperatures as the 
pore size is reduced. 

The capacitance C and tangent loss tan(ô) were measured as a function of 
frequency and temperature for bulk nitrobenzene and for nitrobenzene adsorbed 
in CPG and Vycor glass of different pore sizes ranging from 50 nm to 4.0 nm. 
The typical behavior of x' vs. T is shown in Figure 1b. For pure, bulk nitro- 
benzene, there is a sharp increase in x' at T= 5.6°C, corresponding to the 
melting point of the pure substance. For nitrobenzene confined in CPG, the 
sample is introduced as a suspension of nitrobenzene filled CPG particles in 
pure nitrobenzene, between the capacitor plates. The capacitance measurement 
yielded an effective value of the relative permittivity of the suspension of CPG 
in pure nitrobenzene. Thus k' shows two sudden changes. The increase that 
depends on pore size is attributed to melting in the pores, while the increase at 
5.6°C corresponds to the bulk melting. The shifts in the melting temperature 
obtained using both DSC and dielectric spectroscopy (DS) measurements are in 
good agreement and vary linearly with the reciprocal pore width for large 
pores [5]. The deviations from linearity, and hence from the Gibbs-Thomson 
equation, are appreciable at pore widths as small as 4.0 nm. The DSC scans for 
CCl, melting in porous silica showed similar trends to the melting of nitroben- 
zene (see Figure 5(b)). 


Heat flow, Q/mW 


20.0 -100 oo T/c 100 
(a) (b) 


Figure 1. (a) Representative DSC scans for melting of nitrobenzene in CPG, after sub- 
traction of background signal. Results are for: bulk nitrobenzene, and for CPG samples 
of mean pore size 50 nm, 25 nm and 7.5 nm. Each scan is shifted in vertical scale by 
10 mW for the sake of clarity. (b) Relative permittivity, K' for nitrobenzene in porous 
silica, as a function of temperature for different pore widths 
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The spectrum of the complex permittivity (x',xK'' vs. @ ) is fit to the disper- 
sion relation (equation (3)), to determine the dielectric relaxation time T, 
which gives valuable information about the structure of the condensed phase. 
The frequency range in this study is expected to encompass the resonant fre- 
quencies corresponding to the dielectric relaxation in the solid phases. To probe 
the liquid relaxation behavior would require a frequency range that is 4 to 5 
orders of magnitude higher. The spectrum plots for nitrobenzene in a 7.5 nm 
CPG at temperatures below the freezing temperature in the pore show a Debye 
type relaxation with a single time scale that is estimated to be T = 1.44 ms. At 
temperatures above the pore melting temperature, (e.g., see Figure 2 at T=- 
4°C), the behavior is significantly different. From the double peak structure of 
the K''(@) and the double inflection in the K'(@) curve, two different dielec- 
tric relaxation times are calculated. There is a shorter relaxation time 
T, = 43.6 us, in addition to the longer component T, = 1.7 ms. The longer com- 
ponent relaxation, T, = 1.7 ms, is attributed to the bulk crystalline phase of 
nitrobenzene. The shorter relaxation component, T, = 43.6 us, is attributed to 
the molecular dynamics of the contact layer. The dispersion corresponding to 
the nitrobenzene molecules in the liquid phase occurs in a frequency range that 
is much higher than measured in this study. 


ao 


Figure 2. Spectrum plot for nitrobenzene in a 7.5 nm pore at T = -4°C; this plot yields 
two distinct dielectric absorption regions. The solid and the dotted curves are fits to the 
Debye dispersion relation 


Simulation. The Landau free energy calculations showed that for the case of 
LJ methane and CCl, confined in hard walled and weakly attractive slit pores, 
œ < 085E; , the freezing temperature in the pore was depressed compared to 


the bulk. For most liquids confined in silica based pores the @ values fall in the 
range of the weakly attractive pores [11]. This behavior is consistent with the 
trends observed in the literature and with the Gibbs-Thomson equation. For the 
case of a hard wall pore, the confined system exists as either a liquid or a solid. 
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For a weakly attractive pore that mimics nitrobenzene in silica, the free energy 
surface in Figure 3(a) shows the presence of three phases. Phase A is the liquid 
phase and Phase C is the crystal phase. An intermediate phase B also exists, 
whose structure is plotted in Figure 3(b). The plots represent two-dimensional, 
in-plane pair correlation functions in each of the molecular layers. The pair 
correlation function of the contact layer (the layer adjacent to the pore walls) is 
isotropic, representing a liquid-like layer, while the pair correlation functions 
for the inner layers show a broken translational symmetry corresponding to a 2- 
d crystalline phase; for this system the contact layers freeze at a temperature 
below that of the inner layers. The existence of such a contact layer phase with 
different thermodynamic and structural properties compared to the inner layers 
is consistent with the experimental measurements of the dielectric relaxation 
times for nitrobenzene in CPG and VYCOR. 


80 


60 


A(®) /kT 


(a) (b) 


Figure 3. (a) The Landau free energy for LJ nitrobenzene confined in a model silica 
pore. The three minima correspond to three different phases. (b) The structure of phase 
B, showing that the contact layer is a fluid while the inner layers are frozen 


3.2. Experiments and Simulations on Carbons 

DSC. The evolution in the DSC patterns of different CCl,-ACF samples ob- 
tained during melting and freezing runs has been studied by Watanabe et al. 
[12]. The positions of the peaks in the DSC spectrum were found to be inde- 
pendent of the temperature scanning rate in the range 1.0 to 5 Kmin.”'. The scan 
for bulk CCl, showed three “exothermic” peaks at 242 K, 228 K and 218 K 
corresponding to liquid to metastable fcc-solid phase, fec to rhombohedral 
phase and rhombohedral to monoclinic phase respectively. The observed tran- 
sitions are systematically shifted by 10 K compared to the values found in the 
literature because of supercooling achieved during the freezing run. Freez- 
ing/melting in the confined system occurs at 299 K; the peak positions in each 
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of the DSC scans for melting and freezing runs corresponding to CCl, confined 
in three different samples of ACF (P5 with H=1.1 nm, P10 with H=1.4nm 
and P20 with H = 1.8 nm) showed an upward temperature shift of 57 K. Unlike 
the prediction by the Gibbs-Thomson equation, the freezing temperature is 
found to be independent of pore width in the micropore regime (H = 1.0- 
1.8 nm); see Figure 5. The enthalpies of freezing for the confined system, cal- 
culated from the peak areas were found to be 1.0 % of the bulk values [12]. The 
nature of the solid phase cannot be determined by DSC experiments alone. 
However, the peaks corresponding to the solid-solid transitions that take place 
in the bulk are absent in the confined system. 


Simulation. The Landau free energy calculation for LJ CCl, in ACF's were 
performed for a range of pore sizes. For pore widths that accommodate one to 
three confined CCl, layers, it was found that all the layers froze in unison at the 
freezing temperature, and that this temperature was much higher than the bulk 
freezing temperature of LJ CCl,. For pore widths that accommodate four or 
more layers of adsorbate, it was found that the contact layers froze at a higher 
temperature than the inner layers which in turn froze at a temperature above the 
bulk freezing temperature. The Landau free energy formulation provided the 
means to calculate both transition temperatures. Each Landau free energy curve 
at a particular temperature (see Figure 4) showed a triple-well structure (with 
three minima), corresponding to (1) a phase with isotropic fluid stture within 
each layer, (2) a phase with frozen contact layer (the inner layers being fluid- 
like) and (3) a phase with all the layers having the structure of a 2-d hexagonal 
crystal. 


Figure 4. (a) The Landau free energy for CCl, confined in a model graphite pore. (b) 
The structure of phase B showing that the contact layer is frozen, while the inner layers 
remain fluid-like 
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The cross-over of the grand free energy functions for each phase gave the ther- 
modynamic freezing temperature of the contact layers and the inner layers. The 
freezing temperatures of the contact and the inner layers for the various pore 
widths are shown in Figure 5(a). 


4. DISCUSSION 


The melting point of nitrobenzene and CCl, in the silica pores is always de- 
pressed (Figure 5(b)), due to the rather weak interaction of the fluid molecules 
with the pore walls. The linear relationship between the shift in the pore melt- 
ing temperature and the inverse pore diameter is consistent with the Gibbs- 
Thomson equation for larger pore sizes . The deviations from linearity, and 
hence from the Gibbs-Thomson equation are appreciable at pore widths as 
small as 4.0 nm, however. The quantitative estimates of the rotational relaxa- 
tion times in the fluid and crystal phases of confined nitrobenzene support the 
existence of a contact layer with dynamic and structural properties different 
than the inner layers. The Landau free energy calculation for the simple model 
that mimicked the weak silica wall interaction confirmed the existence of such 
a contact layer with different structural properties. 


Melting of CCI, in Pores 


Freezing of CCI, in ACF: Simulation 


Lay 
425.0 144754 3 2 1 Layers 


375.0 


275.0 


(a) (b) 


Figure 5. (a) Freezing temperatures of the contact layers and the inner layers for LJ CCl, 
confined in ACF slit pores of different pore widths, calculated using the Landau free 
energy method. (b) The shift in the freezing temperature of CCl,, Al = 1, — 1p asa 
function of pore width, in different porous materials obtained from simulation and ex- 
periment 


The importance of the strength of the fluid wall interaction on the nature of 
the shift in the freezing tempreature of fluids in confined systems is shown in 
Figure 5(b). As opposed to the silica based pores, the CCl, molecules experi- 
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ence an elevation in the freezing tempreature when confined in the graphite 
based pores of the activated carbon fibers, consistent with the strong interaction 
potential between the CCl, molecules and graphite. The comparison of the DSC 
results and the simulation (Figure 5(b)) shows that the simplified model of 
spherical LJ CCl, and regular slit shaped graphite pore with smooth walls cap- 
tures the plateau in the A7(H) function; however it overestimates the shift in 
the freezing temperature. In this regime, AT is 60 K in the experiments and is 
about 90 K in the simulation. The overestimation of the freezing temperature 
shift in the simulation is a reflection of the simplicity of the model used. It is 
easier for spherical molecules in slit pores to freeze when compared to five-site 
LJ tetrahedral molecules (a more realistic representation of CCl.) in a real acti- 
vated carbon with a more disordered pore structure. More realistic fluid-fluid 
potential models based on site-site LJ interactions could lead to an improve- 
ment in the prediction of the simulation results, as could a more realistic model 
of the pore structure. The Gibbs-Thomson equation is valid when the effect of 
the contact layers are negligible on the inner layers. When the number of inner 
layers are comparable with the number of contact layers, a deviation from lin- 
ear behavior (G-T regime) is observed. The freezing temperature in the `non- 
linear" regime is influenced by the freezing of the contact layers [9]. 

Since we have calculated the freezing temperatures and the phase behavior 
(Landau free energies) as a function of the relative strengths of the fluid-wall 
and the fluid-fluid interaction, @/€,, a global phase diagram can be con- 
structed that can predict the freezing temperature shift and the phase behavior 
(for example the structure of the different phases in the system) for different 
fluids confined in different porous materials, based on corresponding states 
theory, as long as the interaction energies can be approximated by the LJ po- 
tential [11]. Such a phase diagram can be used as a guide to estimate freezing 
temperatures of different fluids in real porous materials. For example, in the 
case of nitrobenzene confined in ACF, the relative strength of the fluid-wall 
interaction to the fluid-fluid interaction œ is 0.95€, , which corresponds to a 
neutral pore. Thus the corresponding states theory predicts the shift in the 
freezing temperature of this confined system to be close to 0 K [11]. The melt- 
ing temperature of the real nitrobenzene in ACF is about 8°C, measured using 
dielectric relaxation spectroscopy. Thus, the shift is about 2°C, in excellent 
agreement with the predictions of the corresponding states theory [11]. Efforts 
are also underway to use more realistic fluid potentials and pore models in the 
simulation. Recently Gelb and Gubbins [13] proposed a novel mechanism to 
realistically model the formation and structure of porous silica glasses using 
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spinodal decomposition of a binary fluid mixture in the two-phase liquid-liquid 
region. This model is known to closely represent the actual pore size distribu- 
tion of real porous silica and incorporates complex pore networking. We plan to 
study freezing of simple fluids using the free energy method in such a pore 
model, and to make direct comparison with experimental results for CCl, in 
CPG. For the case of activated carbons, Thomson and Gubbins [14] have pro- 
posed an improvement over the slit pore model, by explicitly accounting for 
interactions between different graphene sheets. Their proposed model includes 
networking and polydispersity effects and reproduces the structure of the real 
porous material very closely. 
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grant from the U.S.-Poland Maria Sklodowska-Curie Joint fund (grant no. 
MEN/DOE-97-314). Supercomputer time was provided under a NSF/NRAC 
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